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ABSTRCT 

The focus of the paper is on the derivation of sensitivity equations for transient 
heat transfer problems modeled by different discretization processes. Two examples 
will be used in this study to facilitate the discussion. The first example is a coupled, 
transient heat transfer problem that simulates the press molding process in fabrication 
of composite laminates. These state equations are discretized into standard ^-version 
finite elements and solved by a multiple step, predictor-corrector scheme. The 
sensitivity analysis results based upon the direct and adjoint variable approaches will 
be presented. The second example is a nonlinear transient heat transfer problem 
solved by a p - version time-discontinuous Galerkin’s Method. The resulting matrix 
equation of the state equation is simply in the form of Ax = b , representing a single 
step, time marching scheme. A direct differentiation approach will be used to 
compute the thermal sensitivities of a sample 2D problem. 


INTRODUCTION 

Sensitivity analysis is defined in this paper as a process that derives sensitivity 
equations to compute the derivatives of responses or states with respect to specified 
variables. Since the derivative information can greatly enhance the robustness and 
accuracy of curve fitting, sensitivity analysis becomes a necessary element in many 
engineering applications. Examples include design trade-off, weather prediction, 
analysis error correction, model adjustment, reliability analysis and design 
optimization. 

In sensitivity analysis, the response or dependent variable can be a real number, a 
function or a fimctional and the independent variable can be a real number or a 
function. The challenge of sensitivity analysis arises when the responses to be 
differentiated involve the solutions of some governing equations. In these cases, the 



responses or states are implicitly related to the independent variables through the 
governing equations. These governing state equations can be expressed as 
differential, integral or algebraic forms. The latter is usually as a result of numerical 
discretization of the former. 

Many approaches are available for sensitivity analysis, including automatic 
differentiation 1 , complex variable method 2 , finite differencing, or more traditional 
analytical approaches 3-6 . The analytical approaches used for sensitivity analysis can 
be further classified into various categories. It can be classified as the discrete 
approach or the distributed (continuous) approach, based upon which types of 
governing equations, responses and independent variables are considered. The 
discrete approach works with discretized algebraic equations and real numbers, while 
the distributed approach works with functionals and functions. Sensitivity analysis 
can also be classified as the direct differentiation approach or the adjoint variable 
approach, based upon whether the derivatives of the state variables are computed 
explicitly in the process. The bulk of the effort of the direct differentiation method is 
to establish a sensitivity equation solved for the derivatives of state variables, whereas 
the effort of the adjoint approach is to form and solve an adjoint equation and 
eventually, compute the derivatives of the responses without explicitly computing the 
derivatives of the state variables. 

Subjects related to sensitivity analysis of dynamic transient problems can be 
found in the literature 3 "*' 6 . Thermal transient analysis can be performed in the similar 
fashion 7 " 9 . Recent development and applications of thermal transient analysis can be 
found in aerospace applications 10-11 , laser surface treatment 12 , material processes 1314 , 
and their references. The major focus of the report is on the derivation of sensitivity 
equations for transient heat transfer problems represented in various forms of 
governing equations. Two examples will be used in this study to facilitate the 
discussion. 

The first example is a coupled, transient heat transfer problem that simulates the 
press molding process in fabrication of composite laminates. The state equations are 
made of a heat conduction equation that calculates the through-thickness temperature 
distribution and an empirical equation that monitors the chemical-kinetic reaction of 
resins. These state equations are discretized into standard ft-version finite elements 

and expressed in the form of Ax+ Bx = c . The resultant equations are then solved by 
a multiple step, predictor-corrector scheme. Both of the direct and adjoint variable 
approaches will be used to derive the sensitivity equations in continuous forms. The 
numerical implementation aspects of those sensitivity equations and their accuracies 
will be studied in this paper. 

The second example is a nonlinear transient heat transfer problem solved by a p- 
version, time-discontinuous Galerkin’s Method 15 ' 17 . The resulting matrix equation of 
the state equation is simply in the form of Ax = b, representing a single step, time 
marching scheme. A direct differentiation approach is used to compute the thermal 



sensitivities of a simple 2 D heat conduction problem. Here, possible usages of the 
sensitivity results are presented. 


EXAMPLE 1: A- VERSION FINITE ELEMENT APPROXIMATION 

The interest of this study is in compression molding of a filled polyester resin 
reinforced by chopped glass fibers. The unmolded composite is produced in sheets 
which are from 3 to 6 mm thick, typically. The resin consists of a thick dough and the 
chopped fibers (about 25 mm long) which are randomly oriented in the plane of the 
sheet. In this form, the material is called sheet molding compound, or simply SMC 18 . 
A diffusion reaction system in terms of the temperature distribution and the degree of 
cure can be used to describe the cure process of the composite material under 
consideration. 


The system equations can be expressed as 


du d 2 u da 

pc~ — k— T = pH r —, 
dt dx 2 dt 

in (0,/t)x(0,r] 

(1) 

u{h,t) = u c (t) 

on (0,r] 

(2) 

du(0,t) 

dx 

on (0, t ] 

(3) 

and 



“ = “oM» 

in (0 ,h) at/=0 

(4) 

where p and c are the density and the specific heat of the composite 

material. 


respectively, k is the thermal conductivity in the direction perpendicular to the plane 
of composite material, 2 h is the total thickness, u c (t) is the cure cycle in K ° , T is the 
time in seconds needed for the completion of one cure cycle, H r is the total or 
ultimate heat of reaction, and the last term in Eq. 1, p H r ^ a /^ t , is the rate of heat 
generated by the chemical reaction as characterized by the degree of cure a . 

The degree of cure, a , is defined as the fraction of heat, //(f), released up to time, 
f, for the resin system under cure; a = H(t )/ H r Both H(t) and H r in Eq. 4 can be 
measured experimentally by Differential Scanning Calorimetry (DSC). For an 
uncured material, a approaches zero, and for a completely cured material, a 

approaches one. The reaction rate, . depends strongly on the curing 

temperature. As an example, the cure rate equation of a stepwise isothermal curing 
process which can be used for a polyester SMC is described as follows: 


IT -t<« T > 

dt 

= ( Kj+K 2 a n Xl-a) n (5) 

= ( a 2 e~ d ' ,Kr + a 2 e~ dl /RT a m ){l-af 


where a x ,a 2 ,d l ,d 2 ,m and n are constants, R is the universal gas constant, and K x 
and K 2 are exponential functions of the temperature. 

The optimal cure cycle design 19 aims to select the profile of cure temperature, 
u c (t ) , to achieve the following goals in a compress molding process. 

a) The maximum temperature inside the composite during the cure process can 
not be too high to avoid buring. 

b) The material is cured completely at the end of the cure process. 

c) The material is cured uniformly at any time during the cure process. 

The first two objectives may be mathematically formulated as point-wise functions, 

T(t)<T f in [0,h]x[0,r] 

a(x,T)>a f in [0,/t] 

where x is the total time required to complete one cure cycle. Furthermore, 
the last objective may be measured by the temperature uniformity, which is expected 
to lead to the uniformity of the curing reaction. Mathematically, the temperature 
uniformity is represented by the least-square integral of the deviation between the 
point-wise temperature and the averaged temperature as 


y/ 0 = \[ ju 2 dx-( Judx) 2 / h]dt 
0 0 0 


( 6 ) 


To support the optimal cure cycle design, the thermal design derivatives of the 
functional, y / 0 , the temperature, u(x,t) and the degree of cure, oc(x,t), with respect to 
the cure temperature, u c (t ) , are required. 

Note that the cure temperature appears as a part of the non-homogeneous 
boundary equation in Eq. 2. By using the following replacement of the temperature, 

u(x, t), by u (x, t) as 

u(x,t) = u(x,t)+u c (t) 

which leads to simplification of the heat conduction equation, Eq. 1, 
du , d 2 u du ( - \ 

**= k w- pc ir +pH ' f(a - u+u ' ) 


( 7 ) 


with the homogeneous boundary equations 


u{h,t)~0, 

dx 


in [0,7 ] (8) 

in [0, r ] (9) 


and the initial condition. 


u(x,0) = u 0 (x)-uJ0), in [0,h] (10) 

where /is defined in Eq. 5. Since the initial temperature, u 0 (t ) is the same as the 
initial cure temperature for most applications, Eq. 10 may yield a homogeneous initial 
condition as well. From here on, u(x,t) is abbreviated as u(x,t) for simplification. 


The weak forms of the above equation can now be derived based upon the 
Galerkin’s method for arbitrary functions, w(x) and s(x), as 

hf 3.. ni.. j.. "\ 


n. 


-•-fl 


du d 2 u du 

pc- l c -— + p c £ 

dt Bx 2 dt 


— pHrfydz 


(ID 


and 

^ f)/Y 

n a =0=\(~f)sdx (12) 

o ot 

THERMAL SENSITIVITY ANALYSIS 

The system equation of Eqs. 1-5 simply reveals the fact that their solutions, u(jc,r) 
and a (x,t) are functions of the design variable, u c (t). However, since the design 

variable u c (t) itself is a function, the thermal derivative of u(x,t) with respect to u c (t) 
can be defined as the variation of u(x,t), du , due to the variation in u c (t), 8u c , 

du = ~u( x,t;u c + edu c 

The cure derivative, 6a can be defined similarly. 

The variation of the functional defined in Eq. 6 is then given as 

dyf = \[ j2u - — fudx Jdudxdt 
oo ho 


( 13 ) 



It is assumed that u(x,t;u c ) and a(x,t;u c ) have enough regularity in the time- 
spatial domain and in the design space. Thus, the order of the differentiation and the 
variation is exchangeable; therefore, 

d(&) _ 
dt larj 
d(Sa) _ / 3a N 

dt ^ dt j 
and 


a® 

dx 2 


= J[v* 


With the aid of the above equations, the variations of the state equations defined by 
Eq. 11-12 yield 


T h 


0=11 


00 




A- pH r -^—SaA 
dx 


da 


- pH r f - SuA - pH r 5u c A 

du du , 


dxdt. 


OllU 

0 = ))[sS^^--s^-Sa-s^-du-s^5u\xdt. 
ooV dr da 3u du c J 

Integrating by parts simplify the above equation as 

VJ7 dA ,d 2 A „ ,df\ f dA „ . 


du 


du. 


c / 


-pH r A^-6a 

da 


dxdt+ f k~du dt 
„ J dx „ 


- JfeA d(*^) + dlx-t- Jpc/ldu ( .| dx, 

n dx , n n n n 


(14) 


(15) 


( 16 ) 


and 


( 17 ) 



THE ADJOINT VARIABLE APPROACH 

Adding Eqs. 13, 16, and 17 up, one has the variation of the functional \j/ as 


x Vrl f , 2 \ A 3/1 , d 2 A „ df df 




oo\ 

rh( 




0 0 
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dA 


df 3/ 




Su^dxdt 


cj 


T * 22 * r 2/ii 

+ \pcASu c | + Jjfe— — Su\ dt- \kA— — 
o t o OX 0 0 ox 


dt 


Sudxdt 


h * h r 

+ \pcA5u c \ dx + jsSa\ dx 
o ooo 


(18) 


Note that A(x,t) and s(x,t) are arbitrary functions, and the design derivatives Su and 
Sa are the only two unknowns in the above equation. One may now specify the 
variables X and s in such a way that all of the terms associated with Su and Sa are 
dropped. This can be accomplished by introducing the following adjoint equations for 
X and s as: 


n i^A u 

0 = pc— + k—+pH r 


dA 

dt 


3 / Bf 

du du 




(19) 


and 


n ds .df df 

0= ^: +f>H ^T~ +s T- 
dt da da 


( 20 ) 


with the terminal conditions. 


A( 0 ,t) = 0 , 


in [0, h] 


( 21 ) 


( 22 ) 


s(x,t) = 0, 


and the boundary conditions, 


cU 

dx 


(o,/)=o. 


in [0, h\ 

in [0,7 ] (23) 


/l M = o, 


in [0,7] (24) 


Thus, the combination of Eqs, 18-24 provides a simple formula for the design 
derivative of the functional. 


T h i 


dA 


a/ 








| /Ot/l(o)^u c (o)dx 


(25) 


o 

Equation 25 shows that the design derivative of iff , namely, Sl/f , is a functional 
of state variables a and u, and the adjoint variables A and s. The matrix equations 
which solve the nodal values of a and u can be formed by the standard h - version 
finite element method as. 


Ca+ Ka ~ q 


(26) 


Ma = r 


(27) 


Since the adjoint variables of Eqs. 19 - 20 form an “adjoint” diffusion-reaction 
system similar to the original one, the same numerical scheme used to solved the state 
variables a and u can be extended here to compute the adjoint variables A and s. For 
instance, using the same shape functions of u and a to interpolate the adjoint 
variables X and s obtains the following matrix equation for nodal value of A and s in 
the form of 

CA-KA = q * (28) 


and 

Ms = r* 


(29) 


with proper boundary and terminal conditions. 

In general, the adjoint equation cannot be solved simultaneously alongside with 
the original system equation. Because of the terminal conditions, the adjoint 
equations can be solved by either the backward integration along the real time t - axis 
directly or the forward integration along the artificial time r* - axis, provided that the 


independent variable t is changed to t* ax t* = r-t . However, both approaches 
require the solution of the original system equation known prior to solving the adjoint 
equations. 


THE DIRECT DIFFERENTIATION METHOD 

The direct differentiation method is an approach that differentiates the governing 
equations with respect to a design variable directly. The variation Su in <fy/of Eq. 13 
and 5a can be obtained by solving the equations, 8n a =0 and 8n u = 0, in then- 
weak forms such as Eqs. 14-15, or in their differential forms as 


pc 




dt 


dx 


dt 


da 


+ pH T $-8u + pH r $-8u c 
du du „ 


(30) 


and 


dt da du du„ 


(31) 


where the design variable variation, 8u c , is known. 

Again, Eqs. 14-15 or Eqs. 30-31 can be discretized into finite element matrix 
equations, similar to Eqs. 26-27 as Su and da can be interpolated by the same shape 
functions as u and a . These equations can be symbolically written as 


Cab+Ka b =q b 


(32) 


M a u 


= r b 


(33) 


where 8u and 8a are interpolated by the same shape functions as u and a , in which 
a b and a b are the nodal vectors of Su and 8a . 

NUMERICAL RESULTS 

The initial-value problems of the state variables, their thermal derivatives and 
corresponding adjoint variables are all solved by the computer code, DE 20 . The DE 
program is one of predictor-corrector integration algorithms using Adams family of 


formulas. The truncation error is controlled by varying both the step size and the 
order of the polynomial approximation. The DE program is quite easy to use and has 
the capability to manage moderate stiff equations which happen commonly in the 
problems of chemical kinetics. To maintain a unified accuracy in the analysis, the 
computation of two state variables, namely, the temperature and the degree of cure, 
are subjected to the same error tolerance in this study. 

Note that the coefficient matrices of a b and a h in Eqs. 32-33 are identical to 
those of a and a defined in Eqs. 26-27. Therefore, the same numerical scheme and 
the numerical tolerance can be applied to solve both Eq. 26-27, and 32-33 
simultaneously for state variables, a and a , and design derivatives, a b and a b . In 

this way, the state variables and the design derivatives can enjoy the same numerical 
accuracy. 

Regarding the computational efficiency, it is worthwhile mentioning two notes 
here: 

a) Because the coefficient matrices of Eqs. 32-33 are identical to those of Eqs. 26-27, 
the triangular factorizations of matrices C and M are needed to be done only once. 

b) Compared to the original system equations, the right hand sides of the equations 
for computing a b and a b , such as Eqs. 32-33, may have different frequency contents. 
Thus, to maintain the same numerical accuracy, a small time step At may be required 
for the DE program to solve the pairs (a,a ) and (a b ,a b ), simultaneously. 

Once a and a b are available, the design derivative given in Eq. 13 can be easily 
obtained by numerical integration. Another suggestion is to rewrite the integral form 
of Eq. 13 as a differential equation of Slff given as 


dSy/ 

dt 


= J| 

« o 


Sudz 


(34) 


The above ordinary differential equation of Sy/ b can then be solved simultaneously 
with equations of {a,a ) and (a b ,a b ). In this way, one extra equation of design 
derivative Slff for each design variable is added in the design sensitivity analysis. 
However, the accuracy of Sy/ is secured. Eq. 34 is used to generate the current 
numerical results. 

An example which deals with the cure process of compression molding is 
presented in this section to discuss the numerical accuracy of the adjoint variable 
method and the direct differentiation method for calculating the thermal design 
derivatives. The accuracy of the thermal design sensitivity can be checked by using 
the fundamental definition of design derivatives which can be approximated by the 
finite difference. 


The finite perturbation of the design variable, Ab , is defined as the difference 
between a perturbed design b’ and the nominal design b, i.e., Ab = b* - b. As a 
result, it follows that 

Ay/ = w(b')-w{b) (35) 

= yr Ab 


The above equation provides a simple means to check the accuracy of the design 
sensitivity analysis. Nevertheless, the difficulty of this method is the selection of Ab . 
If Ab is too large, the approximation in Eq. 35 is not valid. On the other hand, if Ab 
is too small, the round-off error in the computation of [y/(b ’ ) - y/{b)\ becomes too 
large to ensure the validity of Eq. 35. 

The first example present here deals with the cure process in which the cure 
temperature of the process is assumed to be a constant temperature. The nominal cure 
temperature is taken as 423° k . According to the approximation defined in Eq. 35, the 
results oi Ay/ and Ab shown in Fig. 1 demonstrate that the thermal design sensitivity 
calculated by the direct differentiation method is more accuracy than the results 
calculated by adjoint variable method. Moreover, by using the direct differentiation 
method, one can also get the time histories of the design derivatives of state variables 
as shown in Figs. 2 and 3. 

EXAMPLE 2: p-VERSION, TIME-DISCONTINIOUS FINITE ELEMENT 
APPROXIMATION 


This example will examine a more general heat conduction equation than the one 
presented in Eq. 1. The governing differential equation is given as 


du 3 3 3 
pc~r~ 2 S — 
dt ■=> j= 1 3fcc, 

u = f(x,t) 


k — 
* ax, 


= Q(x,t), in Qx(0,7°] 
on da u x(0,T] 


3 3 r l 

? X L, —[n,] = q,(x,t), 

'=•/=• dXj 
3 3 0M r i 

X2L,— [n^-hfu-T^), 
1=1 J= l ox, 


on dQ. q x(0,T] 
on 3Q A x(0,r] 


(36) 

(37) 

(38) 

(39) 


and 


u = g(x ) , 


in Q, at t = 0 


(40) 


where the temperature, «(x,t), is the only unknown. Eq. 36 represents an initial- 
boundary value problem with Eqs.37-39 being the temperature, the heat flux 



boundary condition, and the thermal convective condition; respectively, and Eq. 40, 
the initial conditions. It is assumed that the heat source, Q, the prescribed 
temperature, /j the flux, q % the thermal film coefficient, h, and the initial value, g, are 
all with proper regularity. In case of material nonlinearity, the specific heat, pc(u), the 
film coefficient, h(u), and the thermal conductivity, £y(u) are assumed to be functions 
of temperature. 


Let /= 0 on Qu- Otherwise, u in Eqs.36-40 can be replaced by u -f to achieve a 
homogenous boundary condition on As a result, the weak form of Eqs. 36-40 can 
be derived, based upon the Galerkin’s Method for an arbitrary function, w(x), as 


\(pc*^-w+ £ £ ky — ^ )dv+ \huwds 
l dt --U=1 dx t dXj £ 

= jQwdv+ J q s wds+ jhT^wds 
a dQ aa. 


(41) 


In the approach of time-discontinuous Galerkin’s method that is under 
development in this study, the time coordinate is treated as the same as the spatial 
coordinate. The time space is also discretized into elements or intervals. Focusing on 
one time interval, /„ = [t„.j, t „ ), the weak form, Eq. 41, can be extended to the entire 
product domain /„xi2 



n-i a 


du 33 du dw 

— -w+m. — — 

dt 1=1 /=• dx t dXj 


)dv ]dt + J7 jhuwds ]dt 

n - 1 dQ. h 


= j[ jQwdv + jq s wds+ jhT„wds]dt (42) 

n-1 £2 dQ. q 3 Q h 

where, w(x,t) is the testing function. Furthermore, to enforce the continuous 
requirement at the interfaces of time intervals, a weighted integral form of constraint, 
is appended to Eq. 42 

J (/x + u + - pc~u~)w + dv = 0 (43) 

a 


In this example, the temperature, u( x, t), is interpolated by a p- version hierarchical 
basis functions as 


u(x,t) =/ ( x,y,z,t )a 

=(**, ,y)®y(z)®m T a 


( 44 ) 


The symbol, ® , represents the outer tensor product operator and the vectors, yrznd 
0, represent collections of basis functions. Particularly, the through-thickness basis 
functions, yr, are made of the Legendre polynomials of the first kind 15 , the temporal 
basis functions, 0, the integrations of the same Legendre polynomials and the in-plane 
basis functions, chosen for triangular elements, as described by reference 21. 

It is assumed that in this study, the film coefficient h is set to zero and the 
relations between the other material properties, pc and Ky , and the temperature are 

defined in a tabulated form, presented as a result of experiments. Therefore, the 
material properties cannot be explicitly specified as functions of position and time as 
required by integration. An approximation is thus introduced to overcome such a 
difficulty. The standard Lagrange polynomials are used here for this purpose. 

The values of the material properties at the Lagrange points are taken from a 
given material table based upon the values of the temperature found at those points. 
The values of the material properties at elsewhere in an element are then obtained 
through interpolation. In this way, the material properties can be explicitly 
approximated as functions of position and time throughout the problem do main . 

As an example, the material property, say pc, can be interpolated in an element 
(/„ x^,) as 

pc(x,t) = (N c (x, y, z,t)) T £ 

= (N q ,( X ,y)®N cz (.z)®N ct (t)) T £ (45) 

= Z Z Z ( *0, (*, . ( Zj )®N ak (t i )) T S ciJk 

i j k 

where, 5 cijk is a component of the vector, £, which takes the value of the material 
properties found in the material table. 

If the material properties are interpolated linearly through the tabulated data, the 
value of the materials can be represented as, using pc as an example, 

^cijk = a yk + Pijk U ijk 

where a ijk and fi ljk are constants taken from the material table based upon the 
temperature value, u ijk , measured at the location (x t ,y t ,Zj ) and at time t k as 


“yk =zlka 

= (# , y, ) ® WiZj ) ® 6{t k )) T a 



SYSTEM EQUATION 


Once the interpolation functions are selected and the nonlinear material properties 
are approximated as explicit functions of position and time, one can proceed to 
integrate the terms in Eqs.42-43 to construct the equivalent matrices. The resultant 
finite element matrix equation for the time interval, /„ , can then be expressed as 

[c„ («„ )+ K n («») + k )k = 9n + W 'n-1 («„- 1 )]*»_! + (46) 


where the subscript, n, indicates that the associated matrix or vector is evaluated with 
functions defined in time interval, /„ . Note that each term in Eq.46 is corresponding 
to an integral in Eqs. 42-43, which can be evaluated based upon the interpolation 
functions of Eqs. 44-45. As an example, the capacitance matrix is given as 

C. - 

P 

The details of an elemental capacitance matrix, C a ^ is given as Eq. A.l in Appendix. 


The matrix equation derived above is based upon the discontinuous Galerkin’s 
method and will be solved in a time-marching fashion. In other words, the matrix 
equation of Eq. 46 will be solved for a time interval at a time, with a n as unknown 
and a„-j as known quantities. Because its nonlinear nature, Eq. 46 can be solved by 
the Newton-Raphson’s method, which leads to a recursive formula 


[c fl (ar , )+^fc' , )+<k' , )+^k' 1 )+jj < 1 )hi =-k 1 (47) 


and the solution is updated by 



i-l 

n 


+ A< 


(48) 


In the above equation, A a‘ n is the improvement of the solution and R‘ n 1 is the 
residual of the nonlinear equation at the i-l iteration, which is defined as 


R? n-rK 


Moreover, the derivative matrices, J c , Jt and J m are obtained by differentiating 
the coefficient matrices C„ , K„ and M* with respect to the unknown vector a„ , 
respectively. This is usually accomplished at the element level. As an example, the 
derivative matrix, J c , is detailed in Appendix as Eq. A.2. 


SENSITIVITY ANALYSIS 


Since the matrix equation, Eq. 46 , is in the form of a static problem, Ax=b. 
Differentiation of Eq. 46 with respect to a design variable, b, gives a sensitivity 
equation for a nonlinear transient heat transfer problem as, 

(C + K + M + +J C + / t + y„ + ).^- 


fdC dK dM + '] dM~ 
[db db db ) n db 


+ M~ 


fon-i , dq 
db db 


(49) 


The derivative matrices appearing in the right-hand side of Eq. 49 are new for 
sensitivity analysis. Generation of such matrices can be tedious and prone to 
mistakes. A typical derivative matrix, dCldb, is given in Appendix. Fortunately, 
because of the nature of approximation, those derivative matrices can be obtained in 
the same way as the original matrices themselves. This becomes evident by 
comparing Eqs. A1 and A3 in Appendix. Furthermore, it is noted that Eq. 49 is a 
da 

linear equation of — — , which enjoys the same left-hand side coefficient matrices as 
db 

Eq. 46 . Thus, the factored matrices saved from the converged analysis, can be reused 
here to solve the sensitivity equation. That makes the sensitivity analysis 
computationally efficient. 

NUMERICAL RESULTS 

An academic example is used here to verify the derived equations. The space 
domain of the 2-D problem is a 1x1 square, which is discretized into two triangular p- 
version finite elements, as shown in Fig. 4. The heat source term, Q, is specifically 
selected so that the solution of the heat transfer problem, Eq.36, is 

u t ( x , y,t ) = xy(x - l)(y - 1 )t 5 


which gives homogeneous boundary conditions and zero initial condition. Both the 
material properties, pc and isotropic k , are assigned by a linear temperature relation. 


pc,k = 20. + m 



The time step is set to be 1 second and the total operational time, t, is set to be 4 
seconds. 

In the numerical exercise, the orders of in-plane and the temporal polyno mials are 
selected to be 4 and 6, respectively, for temperature interpolation, Eq. 44, whereas the 
orders of in-plane and the temporal polynomials are 3 and 5, respectively, for material 
property interpolation, Eq. 45. The total Lagrange points for the in-plane material 
interpolation is 10 as marked in Fig. 4. Moreover, the error is measured by the L 2 - 
norm 15 as 



dxdy \dt 


1/2 


where e = u-u t and u is the numerical solution. 

At the end of the first time interval, r, = 1 second, the error is in the order of 10" 4 . 
The error is growing with the time. At t 2 = 2 second, the error grows to the order of 
10 1 . At f 4 =4, the error becomes 13.31. The major source of the error is expected 
from deficiency in the in-plane interpolation of materials. The material property is 
linear in temperature. Thus, the exact order of in-plane material interpolation has to 
be 4 which is higher than 3, the order used in the current study. Such error will be 
accumulated from one time interval to the next. The results of temperature 
distribution are shown in Fig. 5. 

Case 1 


The slope of the thermal conductivity-temperature relation is considered as a 
design variable. Since only the matrix K depends upon this special design variable, 
the right-hand side of Eq.49 can be greatly reduced to a single term - (dK / db) a. The 
results of the thermal sensitivity coefficients at times from 1 second to 4 seconds are 
shown in Fig. 6. Note that the thermal sensitivity coefficients is interpolated in the 
same way as the temperature; i.e.. 


du t . , da 

— = Y (x,t ) — 
db A db 


where da / db is obtained from Eq. 49. Comparing with the finite differencing, the 
errors of the thermal sensitivity coefficients are less than 1CT 4 for all the time 
intervals. Note that in this example, thermal sensitivity analysis takes only 18% of the 
time required for one thermal analysis. 


Case 2 


Since in this study, the material property is a distributed function, one may then 
assume that the square slab is made of various materials. In this particular case, the 
value of k at each Lagrange point is determined by it’s own material table. If the slope 
of each material table is considered as a design variable, there are 10 independent 
design variables in total. They are marked in Fig. 4. 

Figures 7 and 8 show the thermal sensitivity coefficients of du/dbj and du /dbg, 
where bj and bg are the slopes of the thermal conductivity-temperature relation at 
Lagrange points 1 and 9, respectively. The figures reveal that the design variable, bj, 
effects the change of temperature along the diagonal line, whereas the design 
variable, b^ does the area off the diagonal line more. Finally, all 10 thermal 
sensitivity coefficients of the temperature at the center point are collected and plotted 
out in Fig. 10. The picture indicates the degree of influence of individual design 
variable on the temperature at the center location at different time. 

CONCLUSIONS 

The paper uses two examples to demonstrate the derivation procedure for thermal 
sensitivity analysis. The continuous approach is used in the first example and the 
discrete approach is used in the second example. It is shown in the first example that 
the direct differentiation method can achieve better accuracy in thermal sensitivity 
analysis than the adjoint variable method. Several authors have similar observation. 
Furthermore, Example 1 shows that the direct results of the direct differentiation 
method, thermal derivatives of the temperature, are very useful in design. This 
particular advantage of the direct differentiation method is also demonstrated in 
Example 2. 

The second example presented here only represented an initial attempt to find the 
thermal sensitivity based upon the p-version time-discontinuous Galerkin’s method. 
Though construction of the matrix equation for thermal analysis is complicated, 
construction of that for thermal sensitivity analysis is rather simple. The resultant 
sensitivity equation is also demonstrated to be computationally efficient. However, 
more works are needed to develop error-control capability for thermal analysis and 
sensitivity analysis to ensure the quality of the p-version time-discontinuous 
Galerkin’s method. 
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Figure 1: Thermal Design Derivatives for Press Molding with Respect to the 

Mold Temperature. 
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Figure 2: Profiles of Thermal Design Derivatives of Temperature with Respect 

to the Mold Temperature. 



Figure 3: Profiles of Thermal Design Derivatives of the State of Cure with 
Respect to the Mold Temperature. 



Figure 4: Meshes and Lagrange Points for Interpolation of Material Property 
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Figure 6: Coefficients of Thermal Sensitivity with Respect to Homogeneous 

Thermal Conductivity; y - . 
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Figure 7: Coefficients of Thermal Design Sensitivity with Respect to Thermal 

Conductivity at Material Lagrange Point 1; _ 
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Figure 8: Coefficients of Thermal Design Sensitivity with Respect to Thermal 

Conductivity at Material Lagrange Point 9; — . 
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Figure 9: Coefficients of Thermal Design Sensitivity with Respect to Thermal 
Conductivity; ; i =1 to 10. 
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